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ABSTRACT 

We study the instantaneous virial balance of clumps and cores (CCs) in three- 
dimensional numerical simulations of driven, magnetohydrodynamical, isother- 
mal molecular clouds (MCs). The models represent a range of magnetic field 
strengths in MCs from subcritical to non-magnetic regimes. We identify CCs at 
different density thresholds in the simulations, and for each object, we calculate 
all the terms that enter the Eulerian form of the virial theorem (EVT). A CC 
is considered gravitationally bound when the gravitational term in the EVT is 
larger than the amount for the system to be virialized, which is more stringent 
than the condition that it be large enough to make the total volume energy nega- 
tive. We also calculate, quantities commonly used in the observations to indicate 
the state of gravitational boundedness of CCs such as the Jeans number Jc, the 
mass-to magnetic flux ratio fic, and the virial parameter amr. Our results sug- 
gest that: a) CCs are dynamical out-of-equilibrium structures, b) The surface 
energies are of the same order than their volume counterparts and thus are very 
important in determining the exact virial balance c) CCs are either in the pro- 
cess of being compressed or dispersed by the velocity field. Yet, not all CCs that 
have a compressive net kinetic energy are gravitationally bound, d) There is no 
one-to-one correspondence between the state of gravitational boundedness of a 
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CC as described by the virial analysis or as implied by the classical indicators. 
In general, in the virial analysis, we observe that only the inner regions of the ob- 
jects (the dense cores selected at high threshold levels) are gravitationally bound, 
whereas J^i <^«r! and fi^ estimates tend to show that they are more gravitation- 
ally bound at the lowest threshold levels and more magnetically supercritical, 
g) We observe, in the non-magnetic simulation, the existence of a bound core 
with structural and dynamical properties that resemble those of the Bok globule 
Barnard 68 (B68). This suggests that B68 like cores can form in a larger MC 
and then be confined by the warm gas of a newly formed, nearby H II region, 
which can heat the gas around the core and confine it. 

Subject headings: ISM : clouds - ISM : globules - ISM : kinematics and dynamics 
- ISM : magnetic fields - turbulence - MHD 

1. INTRODUCTION 

Long standing questions in the process of star formation are the formation mecha- 
nisms, structure, and evolution of molecular clouds (MCs) and their substructure of clumps 
and cores (CCs), which are prone to star formation (e.g., Mouschovias 1995; McKcc 1999; 
Williams et al. 2000; Vazqucz-Scmadcni ct al. 2000; Mac low & Klcssen 2004; Elmegrccn 
& Scalo 2004; Burkert 2006; Ballesteros-Paredes et al. 2007). In MCs, which are observed 
to be turbulent (e.g., Zuckerman & Evans 1974; Zuckerman & Palmer 1974; Larson 1981; 
Blitz 1993; WiUiams et al. 2000; Engargiola et al. 2003; Rosolowsky & Blitz 2005; Koda 
et al. 2006) and magnetized (Heiles et al. 1993; Crutcher 1999; Crutcher et al. 2004), it is 
essential to understand how the fragmentation process occurs. This is the necessary step on 
the way to understanding some of the most important results of the star formation process, 
a) the star formation efficiency, both locally in clouds and globally on galactic scales, b) 
stellar multiplicity and, c) the origin of the Initial Mass Function and its relationship with 
the prestellar cores mass distribution. In the last few decades, several theoretical models 
have been conceived in order to describe the evolving structure of MCs and there has been 
a debate on whether MCs and their substructure of CCs are dynamical out-of equilibrium 
structures or if they evolve in a state of quasi hydrostatic or magnetostatic equilibrium (e.g., 
Mac Low & Klessen 2004; Mouschovias et al. 2006; Ballesteros-Paredes et al. 2007). 

There are two major classes of analytical models used to describe the evolving structure 
of MCs cores, concentrated around non-magnetic and magnetic configurations. In the non- 
magnetic case, Jeans (1902) derived the critical mass (i.e., the Jeans Mass) beyond which 
a system becomes gravitationally unstable against the support provided by the gas thermal 
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motions. This analysis has been expended to include the support provided by isotropic, 
micro-turbulent motions (von Weizsakcr 1943,1951; Chandrasekhar 1951a,b; Sasao 1973). 
Since it was soon recognized that MCs are also turbulent on large scales of the order of 
their own size (e.g., Larson 1981; Basu & Murali 2001; Ossenkopf & Mac Low 2002; Brunt 
2003; Hcycr & Brunt 2004), the assumption of micro-turbulence was replaced by turbulent 
motions that are scale dependent (e.g., Bonazzola et al. 1987; Vazquez-Semadeni & Gazol 
1995). In this class of non-magnetic equilibrium configurations, Bonnor (1956) and Ebert 
(1955,1957) derived analytical solutions for the equilibrium of radial density perturbations in 
a self-gravitating, isothermal ideal gas and which are solutions of the Lane-Emden equation 
(Lane 1870; Emden 1907). Column density observations show that the azimuthally averaged, 
radial density profiles of MCs cores have nearly flat or slow decreasing density in the inner 
parts and a surrounding envelope with a steeper density gradient (Ward-Thompson et al. 
1994,1999; Andre et al. 1996; Bacmann et al. 2000; Shirley et al. 2000,2002; Evans et al. 
2001; Caselli et al. 2002; Young et al. 2003; Kcto ct al. 2004). Bonnor-Ebert (BE) sphere 
models have been applied to radial column density profiles observations of dense (average 
density of ~ 10^ — 10^ cm~^) and cold (T ~ 10 K) cores (Bacmann et al. 2000; Johnstone 
et al. 2000; Alves et al. 2001; Evans et al. 2001; Harvey et al. 2001,2003a,b; Racca et al. 
2002; Tafalla et al. 2004; Lada et al. 2004; Kirk et al. 2005; Kandori et al. 2005) and the 
dynamical evolution of BE spheres has been modeled numerically (McLaughlin & Pudritz 
1996; Ogino et al. 1999; Keto & Field 2005; Aikawa et al. 2005). A study of isothermal 
cores whose structure is a solution to the Lane-Emden equation and which are embedded in 
a larger filamentary structure have been performed by Curry (2000). 

The good agreement in some cases between the models and the observational data has 
lead to the assumption that some cores are indeed BE equilibrium configurations, despite 
the fact that they have, in general, an obvious non-spherical morphology (e.g., Barnard 68; 
it should be noted however that Lombardi & Bertin 2001, showed that the BE instability 
is more related to the core's density contrast than to it's shape). Ballesteros-Paredes et al. 
(2003) showed that azimuthally averaged column densities of dynamical, out-of-equilibrium 
cores formed in a turbulent, isothermal cloud, can be fitted by BE equilibrium profiles. 
Steinacker et al. (2004) argued that projection effects can be misleading about the true 
nature of the density distribution of MCs cores. Recently, Keto & Field (2005) and Keto 
et al. (2006) argued that, not only the density profiles, but also the asymmetric molecular 
spectral lines observations in some cores (e.g., Ba 68, Lada et al. 2003; Redman et al. 
2006) can be explained by non-radial oscillations around an equilibrium BE sphere (see also 
Galli 2005). An intriguing point about BE spheres models is that, as the volume of the 
core decreases and central density and density contrast increases, external pressure should 
increase (e.g.. Figs 1, 2, and 3 in Keto & Field 2005). Yet the derived values of external 
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pressures, precisely from BE modeling of cores, by Kandori et al. (2005) do not show this 
effect (Fig. 7f in their paper). Rather, the external pressure they obtain for a sample of 
cores which are supposed to represent a evolutionary stage of equilibrium BE spheres, seems 
to be independent from the density contrast. 

A second class of models, introduces magnetic fields as an additional supporting mech- 
anism against gravity (e.g., Mouschovias 1976a,b; Mouschovias 1987; Basu & Mouschovias 
1994,1995a,b; Shu et al. 1987; Hujeirat 1998). In the framework of this theory, magnetostati- 
cally supported cores, which can also account for thermal and micro-turbulent support, (e.g., 
Lizano & Shu 1989; McLauglin & Pudritz 1996), increase their masses by accreting neutrals 
that drift across the field lines from the surroundings onto the core's surface. This process is 
called ambipolar diffusion (AD, Mestel & Spitzer 1956). AD will hence tend to increase the 
mass-to magnetic flux ratio in the core until it becomes larger than a certain critical value 
allowing for the denser inner parts of the core to collapse gravitationally. Whereas there is 
no doubt that AD can play an important role provided the ionization fraction inside MCs 
cores is small, the debate has concentrated on whether AD, which acts on it own timcscale 
oitAD ~ 10 tff (McKee 1989; Fiedler & Mouschovias 1992; Ciolek & Basu 2001), where % 
is the free fall time, has time to play any role at all or whether dense cores which form in 
convergent flows in turbulent MCs proceed directly to collapse if their gravitational energy 
dominates all other forms of support (e.g., Balsara et al. 2001; Mac Low & Klessen 2004; 
Mouschovias et al. 2006 and references therein). On the other hand, magnetic field strength 
measurements in cores, which can help constrain the AD timescale and the mass-to mag- 
netic fiux ratio (which is the main parameter in the theory of magnetically mediated star 
formation) of the cores, suffers from a major uncertainty which is due to the cores morpholo- 
gies. Whether the cores are spherical or sheet-like structure can lead to an uncertainty by 
a factor of 1/3 due to their unknown inclinations. Thus, cores that are classified as being 
magnetically supercritical, could be subcritical (Crutchcr et al. 2004). Though observations 
tend to indicate that most cores are supercritical (Crutcher 1999; Sarma et al. 2000; Bourke 
et al. 2001; Crutcher & Troland 2003; Crutcher et al. 2004; Mac Low & Klessen 2004), 
Mouschovias et al. (2006) discuss the possibility that the latter estimates are biased towards 
the very innermost parts of the cores and that the supercritical cores are embedded in sub- 
critical envelopes in agreement with the predictions of the AD theory (Cortes et al. 2005; 
Heiles & Crutcher 2005). 

Whether magnetic or non-magnetic, analytical models suffer the fact that they apply 
only to ideahzed geometries, and neglect the complex interaction of the CCs with their 
environment. These facts, along with some observational and theoretical evidences that 
MCs may have lifetimes not longer than their turbulent crossing time and shorter than 
what would be implied by the AD-theory (e.g., Lee & Myers 1999; Ballesteros-Paredes et al. 
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1999a; Palla & Stabler 1999; Elmegreen 2000; Hartmann 2001; Hartmann et al. 2001) and 
that they must be magnetically supercritical (Nakano 1988), have prompted the development 
of mostly isothermal, magnetized and unmagnetized, turbulent or decaying MCs models in 
order to study the formation, statistics and properties of the dense structures formed within 
them (e.g., Padoan 1995,2001; Mac Low et al. 1998; Ostriker et al. 1999; Klessen et al. 
2000,2001,2005; Smith et al. 2000; Heitsch et al. 2001; Gammie et al. 2003; Bate et al. 
2003; Li et al. 2004; Tilley & Pudritz 2004; Clark & Bonnell 2005; Nakamura & Li 2005). 
The emerging picture from these models is that dense structures (i.e., clumps) form at the 
stagnation points of convergent turbulent flows with a tendency for the densest structures 
(i.e., cores) to form at the intersections of filaments in the clumpy distribution. 

In this picture of the dynamical evolution of MCs, the dynamical properties and detailed 
energy balance of their substructure remains poorly understood. Hunter & Fleck (1982) 
showed that the gravitational boundedness of cores can be influenced by the presence of an 
external velocity field (i.e., a reduction of the Jeans mass). In the context of CCs forming 
and evolving in a turbulent parent MC, a valuable tool to quantify their instantaneous virial 
balance and state of gravitational boundedness is the virial theorem (Chandrasekhar & Fermi 
1953; McCrea 1957; Strittmatter 1963; Zweibel 1990). It has been widely assumed that MCs 
and their substructure are in virial equilibrium (e.g., de Jong et al. 1980; Shu et al. 1987; 
Henriksen 1991). This assumption is also widely spread in observational works (e.g., Larson 
1981; Myers 1983; Solomon et al. 1987; Myers & Goodman 1988a,b; Goodman et al. 1993), 
essentially because their observed kinetic, magnetic and gravitational energies are of the 
same order of magnitude. However, the assumption of the virial balance of CCs in MCs 
have been contested by several authors, particularly for the smallest clouds which arc found 
to be, in general, non self-gravitating (e.g., Carr 1987; Blitz 1987,1993; Maloney 1988,1990; 
Herbertz et al. 1991; Bertoldi & McKce 1992; Hcycr et al. 2001). The observations by 
Swift et al. (2006) of a newly discovered core in L1551 (L1551 MC) suggest that this core is 
highly dynamical and its dynamics is not consistent with a slow gravitational contraction. 
However, it is not completely uncommon to see the detailed virial balance calculated for CCs 
in the observations (e.g., Crutcher 1999; Ward- Thompson et al. 2002;2006). In a study of 
a sample of 27 clouds and CCs in MCs, Crutcher (1999) found that the averaged values of 
the kinetic, magnetic, and gravitational volume energies in the virial theorem verify a static 
form of the virial theorem to 20 %. Individually, the different objects show a large scatter 
and do not satisfy the virial equation by > 30%. Ward-Thompson et al. (2006) calculated 
the detailed virial balance for cores B33-SMM1 and B33-SMM2 in the Horsehead nebulae 
and found them to be out-of-equilibrium and in near virial equilibrium, respectively. Their 
estimates, however, took into account only the volume energies and neglected the surface 
terms in the virial theorem. 
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The aim of this work is to assess, for CCs which form in three-dimensional, isothermal, 
turbulent, magnetized MCs models, whether they are virialized and to quantify the impor- 
tance of the different terms in the virial theorem equation. In § |51 we present the virial 
theorem and discuss the different terms that enter it. The MCs simulations are presented 
in § [3], and in § HJ the clump-finding algorithm used to identify the CCs in the simulations 
is presented. Results of the virial analysis for different MC models are described in § |5l In 
particular, we emphasize on the correspondence between the gravitational boundedness of 
CCs as diagnosed from their virial balance on one hand, on the other by classical gravita- 
tional boundedness indicators such as the Jeans number, the mass-to magnetic flux ratio, 
and the virial parameter, commonly used in theoretical and observational works. In § [6], we 
summarize our results and conclude. 



2. THE EULERIAN VIRIAL THEOREM 

Cast in its Eulerian form (Parker 1979; McKee & Zweibel 1992; McKee 1999), the virial 
theorem (EVT), which is nothing else but a re- writing of the momentum equation after it 
has been dotted by the position vector and integrated over a given volume of interest, has, 
for an object of arbitrary shape, the following expression 

1 •• 1 (i$ 

-Ie = 2 (Eth + Ek- Tth - Tk) + Emag + T^ag + W - - — , (1) 

where Ie = Jy P'^^'^dV is the moment of inertia of the object, r is the distance of each point 
in the object to its center of mass, and V its volume. Eth = | /y PdV is the volume thermal 
energy, with P being the thermal pressure, E^ = ^ Jy rriivfdV the volume kinetic energy, 
Emag = ^ Jy^^dV the volumc magnetic energy, Tth = \ ig^iPhidS the surface thermal 
energy, t^ = \ j^TipViVjhjdS the surface kinetic energy, and Tmag = is^i^ij^jdS the surface 
magnetic energy, where T^j is the Maxwell stress tensor which is given by 

T = — ( B B - i^^I 

4 TT V 2 

where I is the diagonal identity tensor. The effects of the gravitational fleld are described by 
W = — Jy pri{d(j)/dri)dV, where (p is the gravitational potential. This term is not exactly 
equal to the volume gravitational energy because the gravitational potential is a result of 
the distribution of matter inside (e.g., clump or core) and outside the object (e.g., the 
parent cloud), with the essential part of the gravitational acceleration being due to the 
mass contained inside the object. The last term on the right hand side in Eq. [1] is the 
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first time derivative of the flux of moment of inertia through the boundary of the cloud, 
$ = §g pr'^VifiidS. In all the previous quantities, v is the velocity of a point within the 
defined object relative to the velocity of the center of mass of the object. We chose to apply 
the Eulerian version of the virial theorem to the selection of CCs in our simulations simply 
because the numerical grid we are using is Eulerian. Using the Lagrangian form of the 
virial theorem would have forced us to trace fiuid elements, which would have required to 
integrate a passive scalar equation (it would be a method of choice for Lagrangian methods 
such as smooth particle hydrodynamics). Additionally the Eulerian approach, using density- 
threshold defined CCs (i.e., see below), which are Eulerian because they trace a given volume 
rather than a given set of particles, is closer to what is regularly done in the observations. 

Few studies have attempted to evaluate the terms involved in Eq. [T]for CCs formed in a 
turbulent medium (i.e., for non-isolated clouds). Ballesteros-Paredes & Vazquez- Semadeni 
(1997) and Ballesteros-Paredes et al. (1999b) evaluated some of the EVT terms in 2D 
numerical simulations and stressed the importance of the surface term r^. Shadmehri et al. 
(2002) evaluated the different terms of the EVT in the case of dense structures formed in 
3D, MHD simulations with grid resolutions of 100^. These studies made the assumption that 
an object is gravitationally bound if (|W^| > |2 {Eth + Ek - Tk - nh) + Emag + TmagD- More 
recently, Tilley & Pudritz (2004,2006), evaluated the non-time dependent terms in the EVT 
in 3D simulations, and have assumed that the sign of the sum of these terms determines 
whether an object is gravitationally bound or not (i.e., negative is bound). 

Due to the dual role of the surface+volume thermal, magnetic, and kinetic energies and 
also the gravitational field, which can, theoretically, both help confine an object or disperse 
110, we argue that the criterion of instantaneous gravitational boundednes^, is not only 

\W\ > \Qvt\, but also {W + Qvt) < 0, where QvT = 2 {Eth + Ek-n- Tth) + Emag + Tmag- 



^In the case of an idealized axisymmetric filamentary structure, a purely poloidal magnetic field helps 
support the object against gravity while a purely toroidal field tends to aid gravity compress the object. 
The net effect of a helical magnetic field which has toroidal and poloidal components in such configurations 
will be of support/dispersion or compression depending on how much the field is poloidally or toroidally 
dominated (for more discussion on this point see Fiege & Pudritz 2000). 

^In the rest of the text we will not repeat, for commodity, the term 'instantaneous' associated with the 
virial analysis, but it is always assumed to be the case. Also, in the literature it is common to find the term 
'bound' associated with both the compressed/bound case (i.e,. W + Qvt < 0, lower quadrants in Fig. [1]) 
and the virial equilibrium case {W + Qyt = 0, horizontal line in Fig.[T|) case (both cases have a total volume 
energy Efot < 0). In this paper, we make a clear distinction between the two cases by calling the first one, 
indiscriminately, bound or compressed (gravitationally bound, compressed, or captured if the CC is located 
in the lower right quadrant in Fig.[T]) and the second one 'equilibrium' or 'virialized'. Objects located in the 
upper quadrants in Fig. [T] are called indiscriminately 'unbound' or 'dispersed'. 
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If for a given object {W + Qvt) < 0, but \W\ < \Qvt\, the object is instantaneously being 
predominantly bound/compressed by forces other than gravity (i.e., magnetic stresses or 
dynamical compression). Fig. [T] shows, in arbitrary units, a diagram describing the state 
of instantaneous gravitational boundedness of an object according to the balance of the 
different virial terms and which we will apply in the next sections to the identified CCs in 
our simulations. Since the virial theorem ignores the spatial distribution of the different forces 
at work in a CC, its use as a criterion for CCs boundedness gives only correct qualitative 
behavior and scalings that deviate by some factor from the correct quantitative solution 
to the fluid equations (e.g., McCrea 1957; Zweibel 1990). For example, as reminded by 
Zweibel (1990), the use of a constant density model for an isothermal, self-gravitating, non- 
magnetized, and pressure truncated CC in the virial theorem gives results that deviate by 
factors 2-3 from the exact solutions (Spitzer 1978; see also McCrea 1957). Therefore the 
W + QvT = line which in our case draws the boundary between the bound (whether by 
gravity or other forces) and unbound regime may not be a very accurate description of the 
equilibrium of forces for objects which fall very close to the line (i.e. W+Qvt ~ 0). However, 
to this date there is no analytical solution to describe the equilibrium of CCs with complex 
morphologies in the presence of magnetic fields and turbulence (i.e., like the Bonnor-Ebert 
theory for thermal gradient pressure force balancing the gravitational force), particularly 
when the velocity and magnetic fields have complex topologies in and around the objects. 

It is very important to stress that the virial theorem equation represents a snapshot 
of the transformed and spatially integrated balance of the different forces in a CC. The 
virial equation, without any additional physical information such as the local and time 
varying equation of state, rate of turbulent energy injection, and the strength and topology 
of the magnetic field, inside and at the boundaries of the object, can not generally be 
used to make a full prediction about the time evolution of CCs. Such predictions require 
the detailed knowledge on how the (magneto)hydrodynamical quantities evolve over time 
inside and at the boundaries of the objects and the full resolution of the time dependent 
(magneto)hydrodynamical equations. Thus, the diagram in Fig. [1] reflects the instantaneous 
statistics of CCs that are instantaneously bound mainly by gravity, bound mainly by other 
forces, or unbound, and mimics the observational situation where information at different 
epochs can not be obtained. Thus, in the most general case, the usefulness of the virial 
theorem for deriving quantities such as the star formation efficiency (for which a key element 
is the ratio of gravitationally bound objects to non-gravitationally bound objects) does not 
rely on its predictive power for individual objects (some predictions can be made in specific 
situations, e.g., if for example the gas is isothermal, see next paragraphs), but essentially on 
large statistical samples that can be found in each of the quadrants in Fig. [TJ 

Another way of defining the gravitational boundedness of a CC is by requiring that its 
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total volume energy Etot be negative. In simplified systems with no effects of the environment 
(i.e., isolated systems), and for objects that are either in virial equilibrium or where the virial 
balance is tipped in favor of gravity (i.e., W + 2 Eth ^ 0), Etot is guaranteed to be negative. 
An example of such an isolated system is the case of a non-turbulent, non-magnetized gas 
sphere with no external compressive or dispersive agents. For an adiabatic index of the gas in 
the sphere of 7 > 4/3 (e.g., a star) the sphere (i.e., star) is in equilibrium if Eg+2 Eth = (i.e., 
Etot < 0). At equilibrium, the sphere will be located on the intersection of the W + Qvt = 
and |W^|/0yr| = 1 lines in Fig. [T] (in this case Qvt is reduced to the volume thermal energy 
contribution 2 Eth and W to the gravitational energy Eg). The gas may condense in a first 
stage in reaction to a compression (and will have temporarily Eg + 2 Eth < 0, and Etot < 0), 
but as in this case 2 Eth grows faster than the gravitational energy (i.e., thermal gradient 
pressure force grows faster then the gravitational force), the compression will come to a halt 
and the gas will expand {Eg -t- 2 Eth > 0, and still with Etot < 0) and return to equilibrium 
(e.g., Hayashi 1966; Larson 1969; Low & Lynden-Bell 1976; Masunaga & Inutsuka 1999). 
Conversely if the sphere (i.e., star) is perturbed and starts to expand {W + 2 Eth > and 
Etot < 0), \Eg\ decreases slower than 2 Eth and the expansion will be halted and the system 
will contract again and return to equilibrium. In all cases, such a system always has Etot < 0. 

However, if the adiabatic index of the equation of state of the gas is 7 < 4/3 (such 
as for an isothermal, 7 = 1, gas; i.e., like our simulations, and eventually MCs), as soon 
as an object is captured/compressed/bound by the gravitational force (i.e., the lower right 
quadrant in Fig. [T]), the pressure gradient force will not be able to balance the gravitational 
force. In the absence of any other support agent, the CC will be compressed further as 
time proceeds, and if it becomes dense enough, will eventually proceed towards gravitational 
collapse. The same object (i.e., with 7 < 4/3) will expand indefinitely if perturbed by an 
expansion off the equilibrium (even if in an initial stage it still has Etot < 0), as in this 
case \Eg\ decreases faster than 2 Eth- If a CC is found in the lower right quadrant in our 
simulations (i.e., instantaneously gravitationally bound), and eventually in the observation, 
and if the gas is known to be isothermal or nearly-isothermal and if the time varying surface 
terms are assumed not to vary substantially as time evolves, and in the absence of stellar 
feedback, then this CC will very likely remain gravitationally bound until it becomes dense 
enough to undergo gravitational collapse. Thus, the definition of boundedness based on the 
total energy applies only to isolated systems for which the future behavior of the forces can 
be predicted. In non-isolated system, the fact that Etot < and the system be defined as 
'bound' at one instant according to this criterion does not imply that it will remain bound at 
a subsequent epoch. Therefore, in these limiting cases, there is no discrepancy between the 
virial analysis and the Etot argument for objects which are truly instantaneously captured 
by the gravitational force or that are in equilibrium. Both cases can be merged to define 
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the 'term' bound in the Etot argument. However, calling an object that is instantaneously 
being dispersed by the gradient pressure force (i.e., Eg + 2 Eth > 0) 'bound' if it has Etot < 
is only possible if the adiabatic index and eventually its time dependence at each location 
of the object allows such conclusion. More generally, predictions based on the total energy, 
similarly to the virial analysis, require the knowledge of how external pressure and also the 
time varying turbulence and magnetic fields at the boundaries of the object will affect the 
total volume energy. The advantage of using the virial analysis over the total energy is that 
even if it cannot predict the time evolution of the magnetohydrodynamical quantities in 
general systems inside and at the surface of the object, it incorporates an evaluation of the 
instantaneous effects of the surface terms while the total energy criterion does not. 



3. SIMULATIONS 

Aspects of the three-dimensional numerical simulations analyzed in this paper are de- 
scribed in Vazquez- Semadeni et al. (2005a). We use the first 4 simulations summarized in 
Table 1 of that paper, keeping the same nomenclature for the different runs. Here, we recall 
their basic features. The ideal MHD equations are solved using a total variation diminishing 
scheme (TVD), which is a second-ordcr-accurate upwind scheme. Its implementation for 
isothermal flows is described in detail in Kim ct al. (1999). Periodic boundary conditions 
are used in the three directions. The Poisson equation is solved to account for the self-gravity 
of the gas using a standard Fourier algorithm. In order to achieve second-order accuracy 
in time, an update step of the momentum density due to the gravitational force is imple- 
mented, as in Truelove et al. (1997). Following the method described in Stone et al. (1998), 
turbulence is constantly driven in the simulation box and the kinetic energy input rate is 
adjusted as to maintain a constant specified rms sonic Mach number Mg — 10. Kinetic 
energy is injected at large scales, in the wave number range k = 1 — 2. In physical units, 
all 4 simulations have a box size of 4 pc, an average number density of 500 cm~^ (i.e., thus 
the column density is ~ 10^^''^ cm~^ which is similar to that of many MCs except for those 
associated with massive star formation (e.g., Crutchcr 1999), a temperature of 11.4 K, a 
sound speed of 0.2 km s~^, and an initial rms velocity of 2 km s~^ (i.e., therefore the initial 
sonic Mach number is Mg — 10). The Jeans number of the box is Jbox = 4 (i.e., number 
of Jeans masses in the box is Mbox/Mjeans,box — J^ox — ^4). The 4 simulations vary by the 
strength of the initial magnetic field ranging from a magnetically subcritical cloud model 
to a nonmagnetic cloud. The initial magnetic field strength in the box for the subcritical, 
moderately supercritical, and strongly supercritical magnetized models are Bq = 45.8, 14.5, 
and 4.6 fiG, respectively. Correspondingly, the j3 plasma and mass-to-magnetic flux values 
of the box for these runs are, (3 = 0.01, 0.1 and 1, and Hbox = 0.9, 2.8, and 8.8, respectively. 
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Several issues pertaining to the set of simulations presented in this paper such as the 
assumption of isothermality, the neglect of ambipolar diffusion, resolution considerations, 
and the choice of a driven turbulence regime have been discussed in some detail in Vazquez- 
Semadeni et al. (2005a). We complement this discussion by justifying our selection of 
simulations where turbulence is driven on the largest scales. Our motivation for analyzing 
such simulations is due to the increased evidence from recent theoretical work and obser- 
vations that turbulence is most likely driven on large scales in the interstellar medium of 
different galaxies (e.g., Stanimirovic & Lazarian 2001 for the LMC; Dib & Burkert 2004,2005 
for Holmberg II; Koda et al. 2006 in the Galaxy). An analysis of the turbulent velocity struc- 
ture by Ossenkopf & Mac Low (2002) and the filamentary density structure in Taurus by 
Hartmann (2002) also suggest that MCs are driven on scales which are of the order of their 
own size, or larger. 

The nature of the large scale turbulence driver(s) in the ISM remains an unresolved 
issue. A minimum level of turbulence might be due to supernova explosions (e.g., dc Avillez 
2000; Kim et al. 2001; Kim 2004; Dib et al. 2006), but several large scale instabihties are 
suspected to play an additional role in that respect, such as instabilities that occur in spiral 
shocks (Wada & Koda 2004; Bonncll ct al. 2005; Kim & Ostriker 2006; Kim, Kim & Ostriker 
2006; Dobbs et al. 2006), large scale gravitational and thermal instabilities (e.g., Wada et 
al. 2002; Dib & Burkert 2005), generic convergent fiows (Heitsch et al. 2005,2006; Audit & 
Hennebelle 2005; Fohni & Walder 2006; Vazquez- Semadeni et al. 2006a,b) and the magneto- 
rotational instability (e.g., Brandenburg et al. 1995; Selwood & Balbus 1999; Dziourkevitch 
et al. 2004; Dziourkevitch 2005; Piontek & Ostriker 2005). Nevertheless, recent studies also 
show that, in star forming clouds, the feedback from protostellar outfiows (Li & Nakamura 
2006; Matzncr 2007) and H II regions (Krumholz et al. 2006) can, not only enhance the 
level of turbulence inside the clouds, but eventually lead to a self-sustained level of turbulent 
motions. 

The simulations start with a uniform density field, but are evolved for half a sound 
crossing time (the sound crossing time is tg = 20 Myr), equivalent to 5 turbulent crossing 
times (the turbulent crossing time is tc = 2 Myr), before self-gravity is turned on. This 
is a necessary step in order to allow for the full development of the turbulent cascade. In 
some simulations (i.e., those with a weak and a zero magnetic field), CCs form and reach the 
stage of gravitational collapse at epochs which are shorter than a turbulent crossing time. In 
the absence of a prior established turbulent cascade, those CCs could not be considered as 
forming and evolving in a driven medium, but rather in a medium with decaying turbulence, 
which would be contrary to the formation of CCs in a driven turbulence regime we want to 
investigate here. For example, in the subcritical simulation M10J4/3.01, we focus on timestep 
223 which corresponds to the epoch 223 x 0.002 x ts ^ 0.446 ts ^ 4.46 tc after self-gravity 
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has been turned on (i.e., time dumps of the physical quantities are performed at intervals 
of 0.002 ts). For this run the turbulent cascade would have been already established at 
this timcstcp. without any initial driving. However, for the moderately supercritical (i.e., 
run M10J4/j.l), and strongly supercritical (i.e., M10J4/31) and hydrodynamical runs (i.e., 
M10J4/3 inf), we focus the analysis on timesteps 40, 30, and, 20 respectively. This timcsteps 
would correspond to epochs of 0.8, 0.6, and 0.4 tc, respectively. It is clear that at these 
epochs, if no turbulent cascade was already established, the medium in which CCs would be 
forming and evolving in these simulations could not be considered a driven one. 

When applying a virial balance analysis to CCs in simulations such as ours, a crucial 
question is: until which time in the evolution of the simulation does the virial analysis remain 
valid. Due to the limitations of the numerical resolution, the internal dynamics of a collapsing 
object can not be accurately followed if the mass of the object exceeds a Jeans mass and if it 
is not resolved by at least 4 cells in one direction (Truelove 1997). This criterion translates 
into a maximum density {rires ~ 256 n), which when reached locally in one cell, indicates 
that the advection of the flow around this cell is numerically biased. Because we are precisely 
interested in this work in the internal dynamics of the dense structures forming in MCs, we 
choose to present the detailed analysis only for time-steps which precede the appearance of 
any collapsed structure in the simulation box. Fig. 2 in Vazquez-Semadeni et al. (2005a) 
shows the evolution of the maximum density in the sample of simulations we analyze in this 
paper. This figure shows that with a decreasing initial magnetic field strength, CCs tend to 
proceed towards collapse at earlier epochs. It should be noted, however, that the simulations 
are evolved beyond the epoch of appearance of the first collapsed object in the box. The 
global dynamics in the cloud, outside the collapsing CC, remains accurately described by 
the code. Any selection of CCs after the appearance of the first collapsing objects, remains 
valid in order to investigate the statistics of CCs in the box with the caveat that it does not 
properly describe the internal dynamics of the population of cores which are gravitationally 
collapsing. 



4. SIMULATIONS ANALYSIS 

4.1. Clump- finding algoritiim 

We developed a clump-finding algorithm that is based on a density threshold criterion 
and a friend-of-friend approach. The clump finding proceeds as follows: cells that have 
a density higher than the density threshold are selected. Among the selected cells, the 
maximum is sorted out and all cells spatially connected to it are saved in a separate list 
that defines the first clump. The density in the first clump cells are then assigned a value 



- 13 - 



close to zero and the second maximum is searched for in the data cube and the procedure is 
repeated until all cells that define the second clump are found. This is repeated until all cells 
that have densities higher than the selected density threshold are assigned to a clump. Since 
our simulations use periodic boundary conditions, the algorithm also checks for clumps that 
extend across the boundaries of the box. Once the cells of each clump are defined, it becomes 
easy to calculate, for each clump, all the physical quantities that enter the EVT, in addition 
to other physical quantities that characterize its structure and dynamics such as its mass Mc, 
volume Vc, average density ric, turbulent velocity dispersion cr Jf|, and angular momentum Jc- 
The density thresholds rithr used to define the CCs are, in units of the average density (i.e., 
n = 500 cm" -3), 7.5, 15, 30, 60, and 100 rfl. A few time-steps in some runs have additional 
clumps determined with thresholds of 500 n. The latter range of density thresholds (i.e., 
from 3.75 x 10^ to 5 x 10^ cm~'^) mimics a large range of critical densities that are necessary 
to excite various molecular lines in real MCs such as the ^^CO (2-1) and CS (1-0) lines which 
have critical densities of ~ 6.2 and 18 x 10^ cm~^, respectively (Rohlfs &; Wilson 1996). The 
cores in our simulations do reach peak densities of ~ 256 n ~ 1.25 x 10^ cm~^ at the limit 
of the Truelove resolution criterion. This value is close to the critical densities of the N2H"'" 
(1-0, Fi=2-l,F=3-2) and HCO"*" (1-0) lines. However, in pre-collapse situations, very few 
cells that have a density of ~ 250 n are present in the simulation box, and CCs made of cells 
selected at this threshold are entirely dominated by numerical noise. 

In order to lower numerical noise, we transform the surface terms in the EVT into volume 
integrals u=ing the Gau.= theorem (i.e.. f^XdS = /, VXdvfl The volume integral, are 

replaced by summations running over all cells belonging to each object. The time derivatives 
appearing in Eq. [1] are calculated by assuming that the CC, within a time lapse dt, has 
moved over small distances dr = vcm x dt following a uniform linear translation (i.e, with 
no acceleration), where vcm is the velocity of the center of mass of the clump at instant 
t. Since we can not control vcm of each object, we use the smallest time separation of 
time-dumps that is available to us and which amounts to dt = 0.002 tg, where tg = 20 
Myr is the simulation box sound crossing timescale. Thus, the second time derivative of the 



■^The velocity dispersion, CTc, is density weighted, and calculated as being = Jy ni{vi — 
v)'^dxdydz/ Jy nidxdydz, where ni and vi are the local density and local velocity module of each cell, 
respectively, and v is the average velocity module inside the CC. 

"^In the figures shown later, different symbols are used to discriminate between the data points at different 
thresholds. Namely, the 7.5 n threshold level is shown with a (+), the 15 n with a {*), the 30 n with a (o), 
the 60 n level with a (A), and the 100 n level with a (□). 

''We calculate the x-direction component of the gradient of a quantity A at position (i,j,k) as being 
{A{i + k) — A{i — 1, j, fc))/2. The components in the y and z directions are calculated in a similar way. 
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moment of inertia is calculated using estimates of Ie at three different time-steps t — dt, t, 
and t + dt. The flux of moment of inertia d^/dt is an average of the two values calculated 
from the three consecutive time-steps. The typical shifts in one direction are of the order of 
1-2 cells, never exceeding 4 cells. Since the shifts in position associated to shifts from time- 
step t to time-steps t — dt and t + dt generally correspond to non-integer values, the physical 
quantities at the shifted time-steps are obtained by using a trilinear interpolation at the 
new position using all neighboring cells to this new position in the three directions. Hence, 
the approximation we make is that the volume of the cloud Vc, follows a linear uniforno 
motion on distances equal to dr = vcm x dt around the central time-step t. This amounts 
to applying a Galilean transformation to each core, thus preserving the Eulerian nature of 
the virial theorem. 

In order to check the accuracy of our clump-finding algorithm for the computation of the 
different physical quantities, the program was tested against a simplified, albeit unphysical, 
test case problem, consisting of a uniform-density sphere with a pure radial dependence for 
both the velocity and magnetic fields. The values of the density and the amplitudes of the 
velocity and magnetic field are arbitrary. The sphere is placed on a uniform, lower density 
grid with a specified resolution. Fig. [2] shows, in percent, the discrepancy between the 
analytical solution for each of the physical quantities involved in the EVT and the numerical 
solution yielded by our algorithm as a function of the number of grid cells present in the 
diameter of the test sphere. As seen in Fig. [21 32 cells in each direction are needed in order 
to bring the uncertainties below the one percent level. Nevertheless, with more than 4 cells 
per direction, the errors on the EVT quantities are not larger than ~ 15 percent. It should 
be stressed though, that the relative errors shown in Fig. [2] are only order of magnitude 
error estimates on the quantities implied in the EVT, because such uncertainties are model 
dependent. In the case of clumps that have a complex, non-spherical morphology, which is 
the case in our simulations, the uncertainty on the physical quantity will be dominated by the 
uncertainty in the direction with the lowest number of cells. Nevertheless, in order to stay 
consistent with the results of Fig. [2l we must keep in mind that cores with small number of 
cells (typically with less than 3 cells in one direction) can be dominated by numerical noise. 

In our evaluation of the EVT terms, we have neglected the role played by the driv- 
ing force which should also appear on the right-hand side of Eq. [H As the three random 



^This approximation amounts to assuming that the real motion of the object's center of mass: dr = 
1/2 acM dt"^ + Vcm dt (where acM is the center of mass acceleration), over small time intervals dt, is 
dominated by the second term on the right hand side. The acceleration of small, law density CCs is expected 
to be larger than for massive, dense CCs. However, the former ones are those whose quantities are more 
affected by numerical noise. 
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components of the driving force for the whole computational domain have not been stored, 
it is impossible to recover this information since there is no analytical description of this 
force. However, we do not expect this omission to be a problem since in the simulations, 
the driving is performed at the largest scales possible, while the CCs are structures with a 
much smaller size, in each direction, than the box size. Therefore, the driving mainly advects 
the CCs, without greatly impacting their internal dynamics. Nevertheless, we have tested 
this conclusion by running anew a simulation of our sample (the strongly supercritical cloud 
model M10J4/?1, but with different turbulent random seed numbers) and dumped the forcing 
components of the velocity field. From this data, we have calculated the ratio of the kinetic 
energy due to instantaneous forcing {SE'^,) to the bulk total kinetic energy of the cloud E'j^ 
(in E'i^ and the velocity components are bulk velocities: i.e., unsubstracted from the 
center of mass velocity components). For the largest and most massive cloud found in this 
model at t = 30 = 1.2 Myr, this ratio is (9.08,7.07,3.71,3.29,3.35) % at the threshold levels 
of (7.5,15,30,60,100) n, respectively. In general, we find that the percentage of the turbulent 
kinetic energy associated with the forcing in CCs is between 3-9 % of the total kinetic energy 
with an average value around 4%. 



4.2. Calculation of the classical indicators 

For each identified CC, we calculate other quantities commonly used in theoretical and 
observational work to assess the state of gravitational boundedness of the objects, namely, 
the Jeans number Jc, the mass-to-magnetic flux ratio fic, and the virial parameter a^r- It 
should be mentioned at this stage that these indicators, when used in observational studies, 
intrinsically suggest that the thermal, kinetic, and magnetic energies act in support against 
gravity, making two basic simplification to the virial equation: a) The neglect of the surface 
termCl, and b) the assumption that turbulent motions inside the clouds are isotropic. The 
Jeans number is defined as Jc = Rc/Lj^c, where Rc is the characteristic size of the CC, 
Lj,c = i'^cl/GpcY^'^ is the mean Jeans length of the CC and the sound speed. The 
mass-to-magnetic flux ratio {M^/ 4>c) is generally expressed in units of the critical value for 
collapse, calculated, in the linear regime, for a sheet-like structure, (M/0)cr ~ {Ait'^G)~^^'^ 
(Nakano & Nakamura 1978). Krasnopolsky & Gammie (2005) showed that this criterion for 
magnetic criticality holds for turbulent magnetized clouds in the non-linear regime. Note 
that Mouschovias & Spitzer (1976) derived a value of (M/0)cr ~ 0.126 G~^^^ for disks with 



^The concept of magnetic critical mass (similar to the Bonnor-Ebert critical mass in the non-magnetic 
case, includes the surface magnetic term; e.g., McKee et al. (1993). However, the surface magnetic term is 
not commonly estimated in observational studies. 
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only thermal support along field lines, and Tomisaka et al. (1988) found, from extensive 
numerical simulations, that {M/(f))cr ~ 0.12 G"^/^. In principle, the magnetic flux of a CC 
should be computed as 0c = /^^ B n dS, where B is the mean magnetic fleld in the CC and 
n is the unit vector normal to a surface that bisects the CC. In complex geometries such 
as those of the CCs found in our simulations, Sh is difficult to evaluate, particularly for small, 
non-spherical cores which have a limited number of pixels along one direction. Instead, as in 
Vazquez-Semadeni et al. (2005a), we use a simpler approach by defining 0c = vri^^B, where 
Rc is the characteristic CC size. We make the common observational assumption that CCs 
are spherical and compute Rc as being (3 \4/4 vr)^/'^. Alternatively, one could consider that 
Rc = V^^^, or by making the assumption that CCs are fiattened, Rc would be the maximum 
separation between the position of the center of mass and any point (i.e., cell) in the CC. 
Fig. \5\ shows that not all condensations are fiattened and even if some are more fiattened 
than others (e.g.. Fig. [T^ . they are not completely cylindrical objects. The interpretation 
of Jc and fic is that Jc measures if a core is gravitationally unstable with respect to the 
thermal support (i.e., Jc > 1), whereas fic measures the importance of the magnetic support 
against gravity. A core is collapsing when Jc > 1 and /ic > 1, while cores with Jc > 1 and 
/ic < 1 are gravitationally bound but remain in a stable magnetostatic state (under ideal 
MHD conditions such as the simulations analyzed in this work). Cores with Jc < 1 are 
Jeans stable, regardless of the value of fic and are likely to re-disperse due to their internal 
dynamics or loose their identity in subsequent local compressions and dispersions of the local 
medium by large scale fiows. The virial parameter (also called gravitational parameter) is 
calculated as 



""^ " G Mc ' 

where Mc is the mass of the CC and cr is the one-dimensional velocity dispersion inside the 
object, calculated as 



^' = f + (4) 

where cTc is the density weighted turbulent three-dimensional velocity dispersion. The virial 
parameter is often used in observational and theoretical studies in order to measure the 
balance between the CCs self-gravity and their kinetic+thermal energies (e.g., Leung et al. 
1982; Magnani et al. 1985; Keto & Myers 1986; Herbertz et al. 1991; Bertoldi & McKee 1992; 
Falgarone et al. 1992; Dobashi et al. 1996; Yonekura et al. 1997; Kawamura et al. 1998; 
Tachihara et al. 2000,2002; Heyer et al. 2001; Krumholz & McKee 2005). A CC is assumed to 
be gravitationally bound if a^r ^ 1? unbound otherwise. We make the common observational 



-17- 



assumption that turbulent motions are isotropic in the CCs and calculate, in Eq. HI the 
one-dimensional turbulent velocity dispersion as being a^j'i^l'^. An a^j^-Mass relation is 
commonly used in observational studies in order to assess the gravitational boundedness of 
CCs in MCs (e.g., Bertoldi & Mckee 1992; Williams et al. 2000). The a^ir - relation 
is observed to follow a^r oc M^. Tab. [3] displays the values of the exponent e found by 
Bertoldi & McKee (1992) and Williams et al. (2000) for CCs in different MCs. A value 
of e = 2/3 is expected for pressure bounded objects (Bertoldi & McKee 1992). Associated 
to the virial parameter is the concept of virial mass, also commonly used in observational 
studies in order to obtain an estimate of the real mass (e.g., Kramer & Winnewisser 1991; 
Jijina & Myers 1999; Caselli et al. 2002). The virial mass is defined as M^j^ = tt«r ^c- CCs 
whose measured mass is comparable to the virial mass are usually assumed to be in virial 
equilibrium (e.g., Caselh et al. 2002). 

An unavoidable aspect of our simulations is that turbulence inside CCs will be affected 
by the effects of numerical viscosity, particularly for the smallest objects (in addition to 
the discretization effects), which for our numerical code, start to be important at around 8 
grid cells. As an example. Fig. [3] shows the kinetic energy power spectrum for the strongly 
supercritical run at timestep t = 30 = 1.2 Myr. As can be seen from this figure, the turnover 
from the linear regime to the diffusion dominated regime occurs at around wave number 
~ 30 which corresponds to ~ 8 grid cells (i.e., in physical units to sizes of 0.125 — 0.140 
pc). However, the kinetic energy power spectrum does not single out dense structures on 
a given scale. Therefore, it is also important to verify how does the three-dimensional non 
density- weighted velocity dispersion, Onw-, (which also includes a thermal component and is 
only a fraction of percent to a few percents smaller than cxc in all cases) depends on the 
CCs sizes. Fig. H] shows this relatioii^ for the supercritical run at timestep t = 30 = 1.2 
Myr. Despite the discreteness of the data, it is possible to observe a change of regime at 
around Logio{Rc) ~ — 1 which corresponds to ~ 6.4 grid cells. It is difficult to evaluate 
precisely how dissipation will affect the velocity dispersion inside individual CCs because 
of their non-idealized structures. However, from Fig. [3], for scales of ~ 5 cells {k ~ 50), 
the deviation at this location to the inertial range regime is about a factor ~ 3 which is 
roughly the factor by which estimates of the virial parameter will be affected. Note that 
this argument applies essentially to small and dense cores (the dense inner parts of larger 
clumps). For small but low density objects, we find that their virial parameter are much 
larger than unity and thus they are unbound. Applying any correction to take into account 
for the effects of the numerical dissipation of turbulence inside the latter cores will make 



^Thc exponent /3 in the relation (t„^ = i?f is found to he /3 — 0.43 ± 0.06, in very good agreement with 
Larson's (1981) observational result. 
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them even more unbound. 



5. 



VIRIAL BALANCE OF CLUMPS AND CORES 



In this section, we present the results concerning the apphcation of the EVT to our set 
of simulations using the algorithm described in § 14.11 Tab. [T] summarizes the ensemble of 
CCs found in each simulation, for different density thresholds and at different selected time- 



stepqj. At a given time-step, each CC is assigned a number corresponding to its place in the 
cumulative sum of CCs at all thresholds (i.e., for example in run M10J4/5.01 at t = 223, 23 
cores have been found, counting all thresholds. Core number 20 corresponds to the first core 
found at threshold rithr = 30 n). Tab. [2] summarizes the numbering of CCs for time-steps 
which are analyzed in detail in the next sections. 



In the subcritical run M10J4/?.01 {fibox = 0.9 and initial magnetic field strength of 45.8 
/xG), CCs are observed to be very transient (see the animation in Vazquez-Semadeni et al. 
2005a), surviving only for time-spans of 0.01 — 0.075 tg ~ 0.2 — 1.5 Myr before being dispersed 
by their internal dynamics, by large scale flows, or by merging with other clumps, and thus, 
lose their identity. We choose to evaluate the EVT terms in the subcritical run at time- 
step 223 corresponding to t = 8.85 Myr after self-gravity has been turned on, which is the 
epoch at which the most massive clump appears in this run. Fig. O displays two-dimensional 
cuts at the position of the peak density of the most massive object found in the simulation, 
viewed along the three directions of the computational box. Overlayed on the density map 
are the projected velocity (right column) and magnetic (left column) field components and 
density contours which show the extent of the clump when selected at the different density 
thresholds. 



^It is important to note that, for a given set of parameters of the simulations (i.e., same scale of tm'bulence 
driving, initial magnetic field strength, number of Jeans masses and Macli number), the numbers of CCs 
and the mass they contain may slightly change from simulation to simulation due to the different seed of 
the random turbulence driver. An example of this effect is given in Vazquez-Semadeni et al. (2005b). 
Nevertheless, the general properties of the simulations are conserved: That is, no collapsed objects form in 
magnetically subcritical clouds, near critical or moderately supercritical runs will always form a moderate 
number of CCs with a few collapsing ones, and the supercritical and non-magnetic runs will always form a 
large number of CCs with many of them proceeding to gravitational collapse. 




5.1. 



subcritical cloud 
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Fig. [6] displays the virial balance for all CCs selected at this time-step at all thresholds, 
i.e., (PlE/dt^ vs. the right hand side (RHS) in Eq. [TJ Since Eq. [1] is nothing else but a 
re-writing of the momentum equation, all points should, theoretically, fall on the equality 
line. Deviations, which are at most of the order of a factor ~ 5 are due to two effects: 1) 
Discretization errors, due to the limited numbers of cells in the CCs, and 2) the neglect of 
the contribution from the driving force. A legitimate question is: can the small values of 
(PlE/dt^ for small CCs be considered an indication of dynamical equilibrium? An alternative 
view to the fact that (PlE/dt"^ and the RHS are of the same order and which answers the 
latter question is given in Fig. [3 This figure shows that, despite the fact that the moment of 
inertia of the different CCs spans over ~ 10 orders of magnitude, the temporal rate of change, 
normalized to its own value, in percent, i.e., \{dlEl dt) / Ie\, changes by less than one order of 
magnitude, irrespective of their masses or sizes. In fact, the normalized rate of change of the 
moment of inertia of the smallest clumps is observed to he larger, thus indicating their more 
transient nature. In Fig. [HI we evaluate the relative importance of the time derivative of the 
flux of moment of inertia, i.e. 1/2 (d^/dt) of all CCs versus the other terms in the RHS in 
Eq. [H This term is observed to be the dominant one for almost all CCs, particularly for 
the largest ones, and is of the same order of magnitude as d'^lE/dt^ (i.e.. Fig. [H]). This fact 
implies that (dlE/dt) ~ — $, which in turn implies that the time variation of the moment 
of inertia of a clump is essentially determined by the flux through its boundary, which is 
another indication of the dynamical and morphing nature of CCs formed in the simulations. 

Since the time dependent term in the RHS is a purely geometrical term which takes into 
account the morphing nature of the CC and its changing boundary and which adapts such 
that the EVT is always verified, the evolution of a CC will be determined by the non-time 
dependent terms in the RHS of the EVT. Namely, by the competition between confining 
forces (negative energies) on the one hand and dispersive ones on the other hand (positive 
energies). Figs.[TD]a,b,c, and d compare the importance of the individual surface energy terms 
to their volume counterparts, and the total surface and volume energy terms, respectively. 
In general, the surface and volume energy terms are of the same order of magnitude, yet 
with scatter around the equality line, particularly for the case of the larger clumps. The 
difference is more important for the kinetic and magnetic energy pairs than for the thermal 
one. The origin of the differences in the E — t relations is twofold. On the one hand, the 
scatter among different objects is due to the difference of the thermal pressure, velocity and 
magnetic field between the surface of the clumps and their interior. Second, for the same 
object, the difference in the scatter at the different density thresholds is simply the result of 
the energy profile of that object and the position of the surface with respect to the density 
peak imposed by the density threshold. The larger differences observed in the kinetic and 
magnetic energy pairs are due to the anisotropic nature of the kinetic and magnetic surface 
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energy terms and are stronger wherever the velocity or magnetic field do not cross smoothly 
across the surface of the CC. For example, In Fig. [5l it is possible to observe how the velocity 
field penetrates the clump's boundaries at the thresholds of 7.5 n and 15 n becoming more 
parallel to the clump surface at the higher thresholds of 30 n, 60 n and 100 n. 

Fig. [m describes the state of gravitational boundedness of the CCs following the ap- 
proach given in §|2l Objects 1,12,20,22, and 23 correspond to the same condensation seen at 
the thresholds of 7.5,15,30,60 and 100 n, respectively (see Tab 12]). At high density thresholds 
(30,60 and 100 n), the clump has an ellipsoidal structure, with axis ratios of ~ 1 : 2. At 
the lower thresholds of 15 and 7.5 n , the clump is linked by a small bridge to a protu- 
berance to the south which causes its morphology to strongly deviate from an ellipsoidal 
configuration (the bridge is only on one side and not fully seen in the cut at the position 
of the peak density). The outer envelope is being dispersed by forces other than gravity 
{Wi + 6yr,i > and < |6yT,i|) causing the clump to be 'peeled' from the outside. At 
the intermediate threshold levels of 15, 30, and 60 n, the clump is being confined by forces 
other than gravity {W + Qyr < and \W\ < \Qvt)- Material is being redistributed from 
the main clump towards the southern condensation; see, in Fig. [5ti, the projected velocity 
field pointing southward in the bridge. The inner part of the condensation observed at the 
highest threshold level of 100 n (i.e., core 23) is seen to be slightly gravitationally bound. 
However, core 23 has only 141 cells and due to its elongation, is dominated by numerical 
noise along its minor axis. Therefore, the condensation is essentially not bound by gravity. 
This is why this object, the most massive that forms in this simulation, does not proceed 
towards gravitational collapse. The condensation corresponding to clumps (3,14), defined 
at the thresholds of 7.5 and 15 n, respectively, is found to be bound by forces other than 
gravity at the lowest threshold level. This object is a low density clump (i.e., its peak den- 
sity is 29 n) and is very fiattened. Fig. [12] shows density cuts at the position of the clump 
peak density in the three directions of the box along with the projected velocity field. The 
velocity field is compressing the object along its smallest dimensions and would it not be 
for the magnetic support (discussed below), it would end up collapsing into a thin sheet, 
eventually following the description given in Burkert & Hartmann (2004). 

Figs. [12] and [TJ] display the values of the Jeans number J^-, the mass-to magnetic fiux 
ratio (normalized to the critical value for collapse) /ic, and the virial parameter a^r for the 
ensemble of CCs in this simulation (at t = 223 = 8.85 Myr). The Jc and values tend 
to indicate that the most massive condensation (i.e., clumps 1,12,20,22,23) is gravitation- 
ally bound, practically at all threshold levels. The object corresponding to clumps (3,14) is 
marginally bound (Jc < 1 and a^r > 1). However, both objects are observed to be mag- 
netically supported {fic < 1) at all threshold levels, which is in agreement with the virial 
analysis of not being gravitationally bound. It is important to note that the ayir estimates 
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(Fig. [T^ tend to indicate a larger number of gravitationally bound objects with respect to 
their true virial balance or estimates. This is the case for the second most massive object 
in the simulation box (corresponding to clumps 2,13,21 at threshold levels of 7.5,15, and 30 
n) as well as the object corresponding to clumps (4,15 at threshold levels of 7.5 and 15 n, 
respectively). Fig. [11] indicates that the first object is unbound, whereas the second one is 
confined by forces other than gravity. These results tend to suggest that the assessment of 
the dominance of self-gravity in a clump or core using virial parameter estimates is at the 
least dubious and should be used in observational studies with caution. A fit to the a^r — 
data in Fig. [H] yields a relation amr oc , which is, within error bars, in good 

agreement with the a^r — relations found by Bertoldi & McKee (1992) and Williams et 
al. (2000) for CCs in several Giant Molecular Clouds and especially with CCs in Cepheus 
0B3. 

Fig. [15] compares the virial mass M^r to the true mass Mc. The range of masses 
over which M^r ~ Mc is ~ [0.3 — 4] M© and overlaps the one observed in Caselli et al. 
(2002) with a similar scatter around the equality line. We also note a tendency, both in 
the observations and in our results, for the most massive objects [M^ir > 4 M0) to have 

< Mc (i.e., a^ir < 1). Finally we observe a population of small, low density objects for 
which Myir » Mc which indicates their highly dynamical nature. Unfortunately, neither 
in the N2H+ cores sample of Caselli (2002) nor in the ^^CO and C^^O of Bertoldi & McKee 
(1992) and Tachihara et al. (2000,2002) does the CCs masses extend to this low values 
that could enable a comparison with our results. The absence of smaller mass clumps 
in low density tracers observations is very likely due to projection effects. This has been 
demonstrated by Gammie et al. (2003) who compared the mass spectra of three-dimensional 
clumps to the mass spectra of two-dimensional clumps obtained from the projected column 
density map. 

A similarity between the real mass and the virial mass does not imply that the objects 
found in this mass range are in virial equilibrium, but they merely indicate, on the one 
hand that our simulations correctly represent real MCs, and that there is an approximate 
equipartition between the volume kinetic-|-thermal and gravitational energies in that mass 
range though the objects are dynamically evolving and, finally, that if the assumptions (i.e., 
see § 14.21) made in the observations are applied to the CCs found in the simulations, it is 
possible to recover similar results. 
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5.2. Moderate supercritical cloud 

We now turn to the moderately supercritical cloud (i.e., run M10J4/5.1, /Xbox = 2.8 and 
intial magnetic field strength of 14.5 /iG). In this simulation, we observe the formation of 
collapsed objects with high peak densities (ripeafc ^ 5000 n) at three distinctive time-steps, 
i.e., frames 44, 130 and 210, corresponding to the epochs of 1.76, 5.2 and 8.4 Myr after the 
simulation has been started (see animation in Vazquez-Semadeni et al. 2005a). We have 
analyzed the virial balance of CCs in this simulations at frames 30, 40, 50 and 210 (see 
Tab. [1]). Only at time-steps 30 and 40, which precede gravitational collapse in any core, is 
the internal dynamics accurately described for all CCs according to the Truelove criterion. 
The diminishing number of clumps detected at the lower thresholds of 7.5 and 15 n in Tab. [T] 
between frames 30 and 210 clearly shows that clumps are merging to form denser structures. 
An example is shown in Fig. [16] which displays, at t = 40 = 1.6 Myr, the ongoing merger 
of the two most massive cores in the simulation box at this epoch. The merger of these two 
condensations leads to the formation of a massive core which proceeds towards gravitational 
collapse at ~ t = 44 = 1.76 Myr. 

The detailed virial analysis is similar to the one presented for the subcritical cloud. 
The general conclusions drawn from the earlier case such as the verification of Eq. [1], the 
importance of the time dependent and surface energy terms in the virial equation, and 
the more transient nature of the smaller and less dense clumps remain valid for the CCs 
found in this run. However, it is interesting to investigate the energy balance of the CCs at 
the onset of the cloud merger observed in Fig. [161 The energy balance of the ensemble of 
CCs at this epoch is displayed in Fig. [T71 The most massive condensation corresponds to 
clumps 1,17,28,35 at the threshold levels of 7.5,15,30,60 n, respectively, and 38 and 39 at the 
threshold of 100 n which are the two merging cores observed in Fig. [THl The condensation as 
a whole is observed to be unbound at the lowest thresholds (i.e., clumps 1 and 17) but is in 
the process of being compressed by the velocity field (i.e.. Fig. [18]). At higher threshold levels, 
the condensation is gravitationally bound (clumps 28 and 39) or marginally bound (clumps 
35 and 38). On the other hand, the calculated Jc, /ic, and a^r values for this condensation 
(i.e.. Figs. [19] and [20]) suggest that it is gravitationally bound at almost all thresholds (except 
for core 39 which has Jc,39 ^ 1), and close to magnetic criticality (/ic ^ 1)0- 

In contrast to the subcritical run, a larger number of condensations are observed to be 
gravitationally bound (i.e., in the lower right quadrant). The condensation corresponding to 



^°If the characteristic radius of the clumps were calculated as the cubic root of the volume and not as the 
radius of a sphere of equivalent volume, the magnetic flux would be larger by a factor of 3 and as a result 
the clumps would be slightly supercritical /ic, 1,17,28, 35, 38 ^ 1- 
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clumps (2,18,29 and 36) is observed to be gravitationally bound at all threshold levels and the 
inner parts of the condensation corresponding to clumps (3,19,30, and 37) are gravitationally 
bound (cores 30 and 37). Condensation (2,18,29,36), however, is numerically under-resolved 
with a total number of cells of 43 (i.e., at the lowest threshold) with a size of two cells in 
one direction and thus its energy balance is unreliable according to the rough uncertainty 
estimates in Fig. [2l Condensation (3,19,30, and 37) is a well resolved object (i.e., 2794 cells 
at the lowest threshold level of 7.5 n). the virial analysis tends to indicate that the inner 
parts of this condensation are more strongly gravitationally bound (i.e., |M^|/|6yT,37| > 
I W^|/|0VT,3o| > |W^|/|0vr,3,i9|)- On the other hand, the classical indicators tend to show 
that the condensation is more gravitationally bound or marginally bound as a whole than 
its inner parts as indicated by the ayir and Jc estimates {amr,3,i9 < ctmr^so < a«r,37; Jc,3,i9 > 
Jc,3o > Jc,3t) and the fact that the inner parts are less magnetically supported than the 
condensation as a whole according to the /Xc estimates (yUc,37 < /^c,3o < /^c,3,i9)- 

In Fig. [iHl we reproduce the data of Fig. [17] but where clumps are cataloged by whether 
they have Ek — Tk > (triangles), which is an indication of a net dispersive effect by the 
velocity field or Ej. — < (diamonds), which is indicative of a net compressive effect 
by the velocity field. From Fig. [TSl it is interesting to note that unbound CCs, not very 
surprisingly, have Ek — Tk > 0, which is an indication that the velocity field is playing a role 
in the dispersion of these objects. It is also important to note that not all of the CCs that 
are being globally compressed by the velocity field {Ek — < 0) are necessarily dominated 
by gravity as is assumed in some studies (e.g.. Field et al. 2006). 

As in the subcritical cloud case, the ayir estimates tend to indicate a larger number of 
gravitationally bound CCs than that implied by the virial and Jeans number estimates. An 
example is the case of the condensation corresponding to clumps 6,12, and 33 (defined at 
the 7.5,15 and 30 n threshold levels, respectively) which is seen to be gravitationally bound 
according to its a^r values (i.e.. Fig. [20]) whereas Fig. [T7] indicates that this condensation is 
bound by forces other than gravity (outer parts, clumps 6, 22) or unbound (inner parts, core 
33). The relationship between the gravitational boundedness estimator and the real mass 
Mc in Fig. 120] is well fitted by a^r oc M-o-60±o.o3 -^j-^ich is also in very good agreement with 
the Oyir - Mc relations of Bertoldi & McKee (1992) and Williams et al. (2000) (see Tab. E)- 
An additional interesting point, both in this model and in the subcritical cloud model, is 
that the ^^CO like cores (i.e., at the lowest threshold level of 7.5 n) can be described as 
gravitationally bound according to their estimates starting from masses of a few solar 
masses (objects 2,3, and 4 in Fig. and objects 3,4, and 6 in Fig. [201 objects 2 and 5 are 
poorly resolved objects), typically [3-6] Mq. Fig. [2T] displays the same general trends in the 
Mc — Myir relation as in the subcritical cloud model with the exception that the mass range 
over which Mc ~ M^r extends downs to M^r ~ Mc ~ 0.25 M0 and that the data points of 
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Caselli et al. (2002) are better encompassed by our sample of CCs found in this model. 

5.3. Strongly supercritical cloud 

In this section, we discuss the gravitational boundedness of CCs that form in the strongly 
supercritical cloud model M10J4/31 (i.e., weak magnetic field strength of 4.6 /iG; fibox = 8.8). 
A visualization of the time evolution of the dense structures formed in this cloud model (see 
animation in Vazquez- Semadeni et al. 2005a) shows the existence of both collapsed, and 
other, non-collapsing, long-lived condensations. The first of these long-lived condensations 
appears at time-step t = 16 = 0.64 Myr and disperses around frame t = 84 = 3.36 Myr. 
At the epoch at which we perform our analysis (i.e., at time-step t = 30 = 1.2 Myr), no 
collapsing object is present in the simulation box. The clump-finding algorithm finds three 
condensations at all thresholds levels (i.e., 7.5,15,30,60 and 100 n), namely the condensations 
corresponding to clumps (1,12,23,32,37), (2,13,24,33,38), and (3,14,25,34,39). Two conden- 
sations are observed at the first four threshold levels (i.e., clumps 4,15,26,35 and 5,16,27,36, 
respectively), and six condensations are observed at lower threshold levels (i.e.. Tab. [1]). 

Fig. [22] shows the virial balance for the ensemble of CCs at the selected timestep. The 
different condensations show a diverse behavior. The central regions of the condensation 
that contains the highest density peak (i.e., object 1,12,23,32,37) are observed to be gravita- 
tionally bound (i.e., cores 32 and 37) whereas when considered with its larger envelope, it is 
gravitationally unbound (i.e., objects 1,12, and 23). Similarly to the cases discussed in the 
previous sections, the values of the classical indicators, Jc, /ic, and amr for this condensation 
suggest that it is bound at all threshold levels, and magnetically more supercritical when 
considered with its extended envelope (i.e.. Figs. [22] and [H]). The second condensation (i.e., 
clumps 2,13,24,33,38) displays a rather complex behavior. The classical indicators suggest 
that it is gravitationally bound (i.e., estimates in Fig. [2^ or marginally bound (i.e., Jc 
estimates in Fig. [23]) . whereas its virial balance tends to indicate that it is confined by forces 
other than gravity (thermal pressure, dynamical compression and the magnetic force). At 
the second highest threshold level (i.e., 60 n), it is observed to be marginally bound (i.e., 
core 33, IH^I/IOvtI ^ 1-5). This is one of the condensation which are observed to be long 
lived. The third condensation (i.e., clumps 3,14,25,34, and 39) is found to be gravitationally 
bound at all threshold levels in the virial analysis and is one of the objects that later pro- 
ceeds to gravitational collapse. Another condensation exhibits the same trend as the first 
two condensations of having a gravitationally bound or marginally bound central region (i.e., 
cores 26 and 35) and an envelope that is unbound (i.e., clumps 4 and 15). 

We note that in this simulation as well, the ayir estimates (Fig. [21]) tend to suggest a 
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larger number of gravitationally bound clumps than what is indicated by the virial analysis. 
The exponent e in the a^r ~ relation shows in this case a shallower value of e ~ —0.49 ± 
0.03 which is smaller than any of the exponents deduced from the observations (see Tab. [3]). 
The M^ir — Mc relation shown in Fig. [2^ exhibits a larger number of objects that depart from 
the unity line and that are gravitationally bound (i.e., Mmr < Mc) down to masses of ~ 0.6 
Mq, which is not observed in the sample of Caselli et al. (2002). The region where there is 
a good reproduction of the observational data is for masses M^j^ ~ > 2.5 Mq. Even in 
that mass regime, the fact that M^j^ is comparable to the real mass is merely an indication 
of equipartition between the volume kinetic+thermal energies and the gravitational term. 
The assumption that these objects are in virial equilibrium is in disagreement with their 
virial analysis. We also note that the M^r estimates of the objects at the lowest thresholds 
are larger than those at higher threshold levels. A fact reported by Tachihara et al. (2000) 
from their ^^CO and C^'^O survey of cores in Ophiuchus. 



5.4. Nonmagnetic cloud 

In the non-magnetic cloud model (i.e., run M10J4/?inf), the clumps are observed to 
evolve quickly towards gravitational collapse (see Fig. 2 in Vazquez-Semadeni et al. 2005a 
and the animation therein). At t = 30 = 1.2 Myr after the start of the simulation, there are 
already 5 independent cores whose average density is higher than 500 n (i.e.. Tab. [T]). We 
should remind that in our simulations, whenever the density in a cell violates the Truelove 
criterion (Truelove 1997) and becomes higher than ~ 256 n, the dynamics around this cell 
becomes numerically unreliable. Therefore, we show the results of the virial analysis at an 
earlier epoch (i.e., time-step t = 20 = 0.8 Myr), at which only one object has collapsed. 
Fig. [2S] displays the virial balance of CCs at this epoch. The most massive object corre- 
sponds to clumps 1,24,47,70,85, and 88 at the threshold levels of 10,15,30,60,100, and 500 
n, respectively. The inner regions of this object are observed to be undergoing gravitational 
collapse (i.e., cores 85 and 88 since e.g., iWgsl » l^iEk^ss + Eth, ss - ^,88 - nh,88)\) while the 
object as a whole is found to be unbound (clumps, 1,24,47 and 70). Two other condensa- 
tions corresponding to clumps 2,25,48,71, and 86 (at threshold levels of 10,15,30,60, and 100, 
respectively) and 14,37,60, and 83 (at threshold levels of 10,15,30, and 60 n, respectively) 
are also observed to have a gravitationally bound core (i.e., cores 86 and 83, respectively) 
and a more dilute, unbound envelope. 

A peculiar condensation is the one corresponding to clumps (14,37,60, and 83, at thresh- 
old levels 10,15,30, and 60 n, respectively), which is found to be gravitationally bound at all 
threshold levels (Fig. [26]). This object (i.e.. Fig. [27]) is small i.e., 64 cells in total, has 3-5 cells 
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in each direction which correspond to a physical size of 0.046-0.078 pc, has sharp boundaries, 
an average number density of ~ 65 n ~ 3.2 x 10^ cm~^ which is nearly-independent of the 
selected density threshold level, a peak number density of 76.4 n ~ 4 x lO'^ cm~^, a mass 
of ~ 1.5 M0 and a non-thermal velocity dispersion of 0.036 km s^^. This implies a one- 
dimensional velocity dispersion of ~ 0.036/3^/^ ~ 0.02 km s~^, which is roughly one tenth 
of the thermal sound speed. Thus, and despite the fact that this core is highly dynamical 
(see the temporal rate of change of its moment of inertia in Fig. [28|l . very interestingly, all 
of its structural and dynamical properties resemble very closely those of Bok globules and 
particularly those of the Bok globule Barnard 68 (e.g., Hotzel et al. 2002a,b; Lada et al. 
2003). Unfortunately, due do the low number of cells present in each direction, it is not 
possible to draw, for this core, the velocity profile in order to check if they present any of the 
observational evidence of spatial asymmetry in the blue-shifted and red-shifted components 
of molecular lines such as in the case of B68 (Redman et al. 2006). 

One noticeable fact about this non-magnetic cloud simulation is that the number of 
condensations that are found to be gravitationally bound or readily collapsing is not par- 
ticularly larger than the number of gravitationally bound objects found in the magnetized 
runs M10J4/3.1 and M10J4/31. However, though the mass range of CCs in all simulations 
is roughly the same, the total number of CCs and particularly the number of small CCs 
is much larger in the non-magnetized cloud than in the magnetized ones (see numbers in 
Tab. [1]). Figs. [29] and [30l which display the Jc and a^r estimates in this run, respectively, 
illustrate this effect; i.e., same mass range of the clumps as in the magnetized runs but with 
larger numbers of small size CCs. Therefore, the presence of a magnetic field in the cloud 
seems to influence the star formation process by essentially reducing the number of formed 
cores within a certain volume rather than by substantially preventing or delaying the collapse 
process in individual clouds once gravity has overcome all other forces. 

Similarly to the cases of the magnetized clouds, the Jc and a^r estimates for the en- 
semble of CCs in this simulation (i.e.. Figs. [291 and [301) suggest that the most massive con- 
densations (condensation 1 corresponding to clumps 1,24,47,70,85 and 88 and condensation 
2 corresponding to clumps 2,25,48,71,86) are gravitationally bound at all threshold levels, 
being also more bound for the cloud as a whole (i.e., when defined at a lower threshold 
level). This is in disagreement with the virial analysis which suggests that only the inner 
regions of these condensations are gravitationally bound. Also, in this non-magnetic case, 
the a^ir estimates would catalog some CCs as being gravitationally bound (e.g., clump 76) or 
marginally bound (e.g., clump 87), in disagreement with the detailed virial analysis, which 
suggests they are not. We find that the a^r — Mc relation for this model (over-plotted to 
the data in Fig. [30]) is amr oc M~°'^^^°'°^. The exponent, e, is in this case close to the value 
obtained for the strongly supercritical model and is generally in less good agreement with 
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the observations (see Tab. [3]). The Myir — Mc relation for this non-magnetic cloud model is 
shown in Fig. [3T] and displays similar characteristics than the supercritical cloud model. The 
different values of the exponent e between the weak magnetic field (strongly supercritical) 
and non magnetic models on the one hand and the stronger magnetic field cases (subcritical 
and near critical) suggests that the presence of the magnetic field in the cores modifies both 
the density profile and the velocity dispersion-radius relations in the cores upon which the 
value of ayir is dependent. Understanding the exponent of the a^r — Mc relation requires a 
careful study of the radial density profiles of the cores drawn along their principal axis and 
is left to a future detailed study. 

After discussing the virial balance in the sample of magnetized and non-magnetized 
simulations, a few additional points are worth mentioning: a) The maximum density that 
can be reached in our simulations inside a core before the dynamics in these cores violates the 
Truelove criterion is 256 n ~ 1.2 x 10^ cm~^ (we recall n = 500 cm^^) which is close to the 
critical densities for the NH3 and N2li^ molecules of ~ 1 — 2 x 10^ cm^^. All of the cores that 
have peak densities of the order of the critical density for exciting the NH3 and N2H"'" lines are 
gravitationally bound and proceed later to gravitational collapse. Examples are core 38 in 
model M10J4/3.1 and core 37 in model M10J4/?1 and the collapsing cores in the non magnetic 
model. The former two cores have peak number densities of 336.7 and 390.5 n and average 
number densities of 142.7 and 164.2 n, respectively. They possess small amr values in the 
range of 0.15 — 0.20 (see Figs. [20]and[24]) and, thus, are bound according to their a^r values 
as well as the virial analysis. Hence our results tend to show that NH3 or N2H''' cores are 
gravitationally bound. However the number of cells in those cores which have densities larger 
than 256 n is small and the estimates of the virial parameter in such cases are likely to be 
dominated by numerical noise; b) We do not observe, in any of the analyzed simulation, CCs 
that are positioned in the upper right quadrant in the {W + Qvt) — (l^l/|0yr|) diagrams. 
This simply means that the extended distribution of mass around the CCs has very little 
influence on the CCs gravitational boundedness and that the gravitational term W is always 
close to the CCs gravitational energy; c) We observe that some CCs, which are dynamical in 
essence, have W + Qvt ~ and |iy|/|6yT| ~ 1 (i.e.. Figs. [13 and [22]). This could lead them 
to be cataloged as being in a state of magneto-static equilibrium or 'virialized', which is not 
supported by their virial analysis as they have non-zero temporal terms. In reality they are 
just merely 'crossing the line' between the instantaneously unbound and bound regimes. 
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6. CONCLUSIONS AND DISCUSSION 

Understanding the gravitational boundedness of clumps and cores (CCs) in molecular 
clouds (MCs) is a key element in better describing the process of fragmentation of a cloud 
which, itself, is the necessary step on the way of better understanding the process of star 
formation in the dense regions of a MC, the origin of stellar multiplicity, and the initial 
stellar mass function (see, e.g., Padoan & Nordlund 2002,2004; Mac Low & Klessen 2004; 
Ballesteros-Paredes et al. 2007; Klein et al. 2006). 

In this paper, we have analyzed the instantaneous virial balance of CCs that form 
in three-dimensional, isothermal, magnetohydrodynamical, driven, and self-gravitating MC 
models. The analyzed simulations vary by the strength of the magnetic field initially present 
in the box ranging from subcritical to non-magnetic regimes. For each simulation, the 
selection of the CCs and the virial analysis is performed at different timesteps. However, in 
order to follow more accurately the internal dynamics inside the analyzed CCs, results are 
preferentially presented at timesteps that precede the appearance of collapsed objects in the 
simulation box. CCs have been identified by a friend-of-friend algorithm as being an ensemble 
of connected cells that have a density higher than a defined density threshold. We have used a 
range of selected thresholds (i.e., from 7.5 to 100 times the average density) that covers a wide 
dynamical range that mimics the usage of various density tracers sensitive to different average 
gas densities in the observations. Once an object has been defined, we calculate for this object 
all the quantities that appear in the virial theorem in its Eulerian form (Eq. [1]), along with 
other quantities such as its mass, volume, average density, internal velocity dispersion and 
three quantities commonly used in observational studies as gravitational binding indicators, 
namely the Jeans number J^, the mass-to magnetic flux ratio (normalized to the critical 
value for collapse) Hc, and the virial parameter Oyir- In the virial theorem formalism, we 
define a CC as being gravitationally bound when the gravitational term is larger than the 
amount for the system to be virialized, which is more stringent than the condition that it 
be large enough to make the total volume energy negative (see § 2). Our results show that : 

a) CCs are not in virial equilibrium, and are dynamical out-of-equilibrium structures 
as indicated by their complex changing geometry over time and the importance of the flux 
of moment of inertia term in the virial equation, which is a quantitative confirmation of the 
highly dynamical, morphing nature of CCs in a MC. 

b) The surface energy terms are of the same order than the volume energy term (by 
up to a factor ~ 5 smaller or larger) and, thus, are very important in determining the virial 
balance in CCs. Particularly important are the anisotropic kinetic and magnetic surface 
energy terms, which are related to the complex topology of the velocity and magnetic fields 
at the boundary of the CC. We observe that the difference between the volume and surface 
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term is larger whenever the velocity or magnetic field does not cross smoothly the surface of 
the object, due for example to the presence of a high density peak. 

c) We observe that CCs can be either in the process of being compressed by the velocity 
field and in that case they have a net kinetic energy E'^. — < or of being dispersed with 
Ek — Tk > 0. Yet, not all CCs that have — Tk < are necessarily gravitationally bound, 
but are simply in the process of assembling from less dense gas by turbulent ram pressure. 

d) Despite their dynamical nature, some CCs are observed to have W + Oyr ~ 
{GvT = 2 {Eth + Ek- Tth - Tk) + Emag + Tmag) and |W^|/|©yr| ~ 1 which can lead them, 
erroneously, to be cataloged as being in a state of virial or magnetostatic equilibrium or 
'virialized', despite their dynamical nature. Additionally, we find, in all simulations, that 
the gravitational term must always be close to the gravitational energy and that the mass 
distribution outside the CCs docs not have a significant contribution in distorting the CCs 
structure by external gravitational torques. 

e) There is no one-to-one correspondence between the state of gravitational boundedness 
of a CC as described by the virial balance analysis (i.e., gravity versus other virial energy 
terms) and as implied by the classical gravitational binding indicators Jc, fJ^c, and a^r- In 
general, from the virial analysis we observe that only the inner regions of the clumps (i.e., 
the dense cores selected at high density thresholds) are gravitationally bound, whereas Jeans 
number estimates of the same clumps tend to show that the objects are gravitationally bound 
at all threshold levels. On the other hand, the calculated a^j^ values not only shows that the 
clumps are more gravitationally bound at the lower threshold levels as the Jeans numbers, 
but also indicate a number of gravitationally bound objects always in excess of what is 
yielded by the virial analysis. Preliminary results by Dib & Kim (2006) show that the Jc 
estimates and the corresponding energy ratio, when including the thermal surface energy 
{Effi — Ttfi/\W\) are well correlated, and that the virial parameter a^r and the corresponding 
energy ratio {Ek + Eth — Tk — Tth)/\W\) or {{Ek — Tk)/\W\) are poorly correlated. A future 
study (Dib et al., in preparation) will investigate these relations in more detail for different 
models as well as the correlation between /ic and the magnetic energy to gravitational term 
ratio. 

f) The exponent, e, of the virial parameter-mass relation, ayir oc M^, shows a better 
agreement with observationally derived values for the case of the moderately supercritical 
cloud in which the most massive CCs are near critical {/ic ~ 1) , and a lesser good agreement 
for the subcritical (in which most CCs are subcritical) and strongly supercritical (in which the 
most massive object is supercritical) and non-magnetic cases. This result is not extremely 
conclusive yet because of the large uncertainties in the observational values of e and the 
limited statistics in our simulations. Yet, it might refiect an agreement with the recent 
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predictions that most observed CCs have a mass-to magnetic flux ratio that is supercritical 
or subcritical by factors of 2 with the median value being around /^c ~ 1 (Crutcher & Troland 
2006). 

g) In the non-magnetic simulation, we observe the formation of a core which possesses 
all the structural and dynamical properties of the Bok globule Barnard 68 (B68). This core 
is gravitationally bound. Such bound cores may survive the ionizing front of a H II region 
formed elsewhere in the cloud and which expands in the cloud's clumpy distribution (e.g., 
Mellema et al. 2006; Will Henncy, private communication). The ionizing front can evacuate 
the gas around the core, and leave it confined by a surrounding warm gas such as in the case 
of B68. 

In the set of simulations presented in this work, we have used an isothermal equation 
of state to describe the gas physics. However, several authors have argued recently that 
deviations from isothcrmality in molecular clouds might lead them to have distinct structural 
and/or dynamical properties than isothermal clouds (e.g., Larson 2005). Li et al. (2003) and 
Jappsen et al. (2004) showed that the fragmentation process is dependent on the polytropic 
exponent which describes the equation of state. Dib et al. (2004) showed that the average 
density-size relation for non-isothermal clouds follows the observed n oc Larson relations 
(Larson, 1981) unlike the isothermal clouds and their substructure is which this relation is 
not foimd (e.g., Ballesteros-Paredes & Mac Low 2002; Li et al. 2004). Hcnnebelle & Inutsuka 
(2006) discuss the possibility that warm gas can survive inside MCs. However, Pavlovski et 
al. (2006), showed that deviations from isothcrmality arc only visible in terms of the high 
temperature zones behind shock waves. The detailed effects of a non-isothermal equation of 
state on the properties of individual cores has not been investigated in detail so far, and is 
left to a future work. 
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Table 1. Number of clumps and cores found in the simulations for the various density 

thresholds at selected time-steps. 



Model 


time-step nth 


r = 7.5 n 


nthr = 15 n 


nthr = 30 n 


nthr = 60 n 


nthr = 100 n 


nthr = 500 n 


M10J4/3.01 


223 


11 


8 


2 


1 


1 





M10J4/J.1 


30 


19 


18 


8 


5 


4 







40 


16 


11 


7 


3 


2 







50 


14 


14 


9 


2 


1 


1 




210 


12 


9 


9 


4 


3 


2 


M10J4/?1 


20 


16 


14 


7 


2 


2 







30 


11 


11 


9 


5 


3 







40 


15 


13 


8 


6 


1 







50 


14 


10 


5 


3 


3 


1 




60 


16 


15 


5 


4 


3 


1 




70 


15 


12 


7 


5 


3 


1 


M10J4/3inf 


20 


23" 


23 


23 


15 


3 


1 




30 


19" 


19 


19 


12 


12 


5 




40 


22" 


22 


22 


18 


12 


6 


Note. — (a) 


due to the large 


number of cells found in 


dense structures 


in the non-maf 


;netic run M10J4/3inf, the lowest 



density threshold used is 10 n. 



Table 2. Numbering of the condensations found in run M10J4/5.01 at t=223=8.92 Myr, run M10J4/5.1 at t=40=1.6 
Myr, and run M10J4/91 at t=30=1.2 Myr. Each column corresponds to a separate object which is attributed a 
distinct number at the different density thresholds. In the case the inner parts of an object are fragmented into 2 
clumps, each core is assigned a different number at the corresponding higher density threshold. The numbering of the 
non-magnetic run is not shown due to the large number of columns needed (i.e., 23 different objects, see Tab. [I]), the 
numbering of CCs in this model follows the same rules as for the magnetized runs. Between parentheses are shown 
the cubic roots of the total number of cells present in each object. This only represents the characteristic number, not 
the number of cells along each of the object's principal axis since many objects are elongated. In our scaling, one cell 

represents a physical scale of 4/256=0.015625 pc. 
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M10.T4^.01 at t=223=8.92 Myr 


























7.5 n 


1 (53.0) 


2 (14.9) 


3 (16.1) 


4 (15.2) 


5 (1.7) 


6 (11.51) 


7 (3.8) 


8 (2.9) 


9 (7.8) 


10 (1.2) 


11 (1) 




15 n 


12 (26.1) 


13 (12.0) 


14 (14.1) 


15 (15.4) 


16 (1.7) 


17 (11.0) 


18 (3.8) 


19 (2.9) 










30 n 


20 (19.0) 


21 (4.8) 






















60 n 


22 (12.4) 
























100 n 


23 (5.2) 
























M10J4/3.1 at t=40=1.6 Myr 


























7.5 n 


1 (58.9) 


2 (3.5) 


3 (14.0) 


4 (16.0) 


5 (6.4) 


6 (18.5) 


7 (7.8) 


8 (10.2) 


9 (11.5) 


10 (1) 


11 (2.1) 12 (1) 13 (3.8) 


14 (2.3) 15 (1.2) 16 (1.4) 


15 n 


17 (33.2) 


18 (3.5) 


19 (14.0) 


20 (16.0) 


21 (6.4) 


22 (18.5) 


23 (7.8) 


24 (10.2) 


25 (7.0) 


26 (1) 


27 (2.1) 




30 n 


28 (19.1) 


29 (3.5) 


30 (10.4) 


31 (9.5) 


32 (6.4) 


33 (9.6) 


34 (7.6) 












60 n 


35 (13.8) 


36 (3.5) 


37 (5.1) 




















100 n 


38,39 (9.8,3.5) 
























M10J4/31 at t = 30=1.2 Myr 


























7.5 n 


1 (59.1) 


2 (12.1) 


3 (19.9) 


4 (11.5) 


5 (16.7) 


6 (7.0) 


7 (1.2) 


8 (2.2) 


9 (11.7) 


10 (5.6) 


11 (1.2) 




15 n 


12 (32.3) 


13 (12.1) 


14 (18.5) 


15 (11.5) 


16 (16.7) 


17 (7.0) 


18 (1.2) 


19 (2.2) 


20 (11.7) 


21 (5.6) 


22 (1.2) 




30 n 


23 (18.3) 


24 (12.1) 


25 (2.7) 


26 (8.4) 


27 (12.4) 


28 (7.0) 


29 (1.2) 


30 (2.2) 


31 (3.3) 








60 n 


32 (10.3) 


33 (10.2) 


34 (7.8) 


35 (2.6) 


36 (3.5) 
















100 n 


37 (7.8) 


38 (4.1) 


39 (3.2) 
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Table 3. Exponent, e, of the Virial parameter-Mass relation, ayir oc obtained from the 
observations of a number of molecular clouds (Bertoldi & McKee 1992"; WiUiams et al. 

2000) and from our models. 



Cloud /Alodel 



G216-2.5 (Williams et al.) ~ -0.72 

Rosette (Williams et al.) —0.61 

Rosette (Bertoldi & McKee) -0.50 ± 0.10 

Ophiucus (Bertoldi & McKee) -0.68 ± 0.03 

Orion B (Bertoldi & McKee) -0.67 ± 0.24 
Cepheus 0B3 (Bertoldi & McKee) -0.54 ± 0.08 

model M10J4/J.01 -0.57 ± 0.03 

model M10J4/3.1 -0.60 ± 0.03 

model M10J4/J1 -0.49 ± 0.03 

model MIO J4/3 inf -0.48 ±0.02 



Note. — (a) Bertoldi & McKcc (1992) did not in- 
clude the massive objects in their fits to the 
relations which would have the effect of giving, for 
some clouds they have analyzed, a slightly shallower 
value of e 
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Fig. 1. — A sketch showing the instantaneous gravitational boundedness state of an 
object according to its virial balance, shown here in arbitrary units. If the cloud in 
which the object forms is turbulent and magnetized, the general form of Qvt is Qvt = 
(2 {Eth + Ek-Tk- Tth) + Ernag + Tmag) , whcrc the Ei and Tj (with i = th, k, mag) terms are 
the volume and surface, thermal, kinetic and magnetic energies, respectively. The detailed 
expression of each term is given in §[2l In the text, we use indiscriminately the terms 'bound' 
and 'compressed' for objects located in the lower quadrants, and the terms 'unbound' and 
'dispersed' for objects located in the upper quadrants. Objects that fall on the horizontal 
W + QvT = are called indiscriminately 'instantaneously in equilibrium' or instantaneously 
virialized'. In isothermal systems such as the cores and clumps in our simulations, the notion 
of instantaneous (un) boundedness can be replaced by a notion of 'tendency to be perma- 
nently bound/unbound'. Thus, in isothermal systems an object located in the lower right 
quadrant will have a tendency of being permanently gravitationally bound by gravity and 
eventually proceed towards gravitational collapse (see §2). 
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Fig. 2. — Error estimates for the terms involved in the Eulerian Virial Theorem for a test 
setup. The test setup is a sphere with a uniform density distribution. The velocity and 
magnetic fields possess, both, a pure (unphysical) radial dependence with arbitrary ampli- 
tudes which facilitates the comparison to the analytical solution. Ek, Eth, E^ag, and Tk, Tth, 
Tmag are the volume and surface kinetic, thermal, and magnetic energies, respectively./ is 
the moment of Inertia and W the gravitational term. 
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0.0 0.5 1.0 1.5 2.0 2.5 

Log^o(^) 

Fig. 3. — Kinetic energy power spectrum for the supercritical run M10J4/31 at i = 30 = 1.2 
Myr. The turnover to the diffusive regime occurs at wave number A; ~ 30 equivalent to ~ 8 
grid cells. 
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-2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 

Logio Rc (pc) 

Fig. 4. — Velocity dispersion-Size relation for the ensemble of CCs in run M10J4/31 at 
t = 30 = 1.2 Myr. The quantity anw includes a thermal component and a three-dimensional, 
non-density weighted, turbulent component. This relation is shown in order to evaluate at 
which scale the effects of numerical diffusion become important for the CCs. 



Fig. 5. — Density cut in the x (top), y (middle), and z (bottom) directions at the position 
of the density maximum for the densest clump in run M10J4/3.01 aX t = 223. This object 
corresponds to clump 1, 12, 20, 22 and 23 in Tab. [H The contours show the boundary of the 
object at the density thresholds of 7.5, 15, 30, 60 and 100 n. The cross marks the position 
of the density maximum. Arrows show the projected velocity (left) and magnetic (right) 
fields at the position of the two-dimensional cut, scaled to the maximum value on the grid. 
In order to better highlight the structure of the cloud, the resolution of the map has been 
artificially multiplied by a factor 2 and the map smoothed. 
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Fig. 6. — Right Hand Side terms in Eq. [T] (RHS) vs. the Eulerian second time derivative 
of the moment of Inertia (flE/df in run M10J4/?.01 at t = 223 = 8.92 Myr. The virial 
theorem is verified to within ~ 0.5 dex for the largest clouds. For the smaller clumps (i.e., 
with small numbers of cells), the larger scatter is due on the one hand to the effects of the 
turbulent driving, unaccounted for in Eq. [Hand to numerical noise (see Fig. [1]). The 7.5 n 
threshold level is shown with a (+), the 15 n with a (*), the 30 n with a (o), the 60 n level 
with a (A), and the 100 n level with a (□). 
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Fig. 7. — This plot shows the temporal rate of change, in percentage, of the moment of 
inertia for the ensemble of clumps and cores in run M10J4/30.01 at t — 223 = 8.92 Myr. 
The smallest clumps exhibit larger temporal variations of their moment of inertia, which is 
indicative of their more dynamical nature. 
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Fig. 8. — Comparison of the temporal term d^/dt to the other terms in the RHS of the 
EVT. The figure is for run M10J4/30.01 at t = 223 = 8.92 Myr. 
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Fig. 9. — Contribution of the d^/dt term to the total second time derivative of the moment 
of inertia. The figure is for run M10J4/30.01 at t = 223. 
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Fig. 10. — Surface energy terms versus volume energy terms a) kinetic, b) thermal, c) 
magnetic and d), total for run M10J4/3.01 at i = 223 = 8.92 Myr. 
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Fig. 11. — This figure compares the importance of the gravitational term (y axis) to the 
net amount of thermal, kinetic, and magnetic energies that appear in the EVT (x axis) in 
run M10J4/30.01 at t = 223 = 8.92 Myr. 




Fig. 12. — Density cuts in the x (top), y (middle), and z (bottom) directions at the position 
of the peak density for the condensation corresponding to clumps number 3 and 14 in run 
M10J4;g.01 at t = 223 = 8.85 Myr, defined at the nth = 7.5 n and 15 n threshold levels, 
respectively. Arrows represent the projected velocity field. In order to better highlight the 
structure of the cloud, the resolution of the map has been artificially multiplied by a factor 
3 and the map smoothed. 
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Fig. 13. — Jeans numbers Jc and mass-to-magnetic flux ratio (normalized to the critical value 
for collapse) of the ensemble of clumps and cores found in run M10J4/3.01 at t = 223 = 8.92 
Myr. 
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Fig. 14. — The virial parameter as a function of the mass of the clumps and cores in run 
M10J4/3.01 at i = 223 = 8.92 Myr. 
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Fig. 15. — The real clumps and cores mass Mc is compared to the virial mass My^^ which 
is calculated by comparing the clumps gravitational and kinetic+thermal energies in run 
M10J4/3.01 at t = 223 = 8.92 Myr. Thick triangles and thick diamonds correspond to the 
starless CCs and CCs which contains stars in the sample of Caselli et al. (2002), respectively. 



Fig. 16. — Density cut in the x-direction showing two merging core cores in run M10J4/?.l 
at t = 40 = 1.6 Myr, respectively. Arrows represent the projected velocity field scaled to its 
maximum value on the map. 
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Fig. 17. — Same as Fig. [TT] but for the moderately supercritical run M10J4/5.1 at t = 40 = 1.6 
Myr. 
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Fig. 18. — Same as Fig. [T7]but where clumps are cataloged by whether they have (E'^-Tfe) > 
(triangles) or {Ek — r^) < (diamonds). A positive/negative (E'^ — r^) is indicative that the 
clump or core is experiencing a net dispersive/compressive effect of the velocity field. 
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Fig. 19. — Jeans numbers Jc and mass-to-magnetic flux ratio (normalized to the critical value 
for collapse) /ic of the ensemble of clumps and cores found in the moderately supercritical 
run M10J4/3.1 at t = 40 = 1.6 Myr. 
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Fig. 21.— Same as Fig. [15] but for run M10J4/5.1 at t = 40 = 1.6 Myr. 
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Fig. 22. — Same as Fig. [TT]but for the strongly supercritical run M10J4/31 at t = 30 = 1.2 
Myr. 



-67- 



.0 



0. 



0.0 



?l 15 

O )K = 

26 



□ 



\1Z 10 36 



0.1 



3* 



.00 



Fig. 23. — Same as Fig. [13] but for the strongly supercritical run M10J4/31 at t = 30 = 1.2 
Myr. 




Fig. 24. — Same as Fig. [TUbut for the strongly supercritical run M10J4/31 at t = 30 = 1.2 
Myr. 
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Fig. 25. — Same as Fig. [15] but for the strongly supercritical run M10J4/31 at t = 30 = 1.2 
Myr. 
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Fig. 26. — Same as Fig. [TTlbut for the non-magnetic run M10J4/5inf at t = 20 = 0.8 Myr. 
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Fig. 27. — Cut at the position of the density maximum in the x-direction of the box for the 
object corresponding to clumps (14,37,60, and 83 at the thresholds levels of 10,15,30 and 60 
n, respectively) in the non-magnetic run M10J4/3inf at t = 20 = 0.8 Myr. The dynamical 
and thermodynamical properties as well as the bean-like morphology of this object resemble 
those of the starless core Barnard 68 (Alves et al. 2001a,b). In order to better highlight the 
structure of the cloud, the resolution of the map has been artificially multiplied by a factor 
8 and the map smoothed. 




Fig. 28. — Same as Fig. [7] but for the non-magnetic run MIO J4/5 inf at t = 20 = 0.8 Myr. 
The Barnard 68-hke object corresponding to clumps (14,37,60, and 83 at the threshold levels 
of 10,15,30 and 60 n), shows a temporal rate of change of it's moment of inertia of the order 
of ~ 10%. 
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Fig. 29. — Jeans number-Mass relations for the ensemble of clumps and cores in run 
M10J4/3inf at i = 20 = 0.8 Myr. 




Fig. 30.— 



Same as Fig. [H but for run M10J4/5inf at t = 20 = 0.8 Myr. 
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Fig. 31.— Same as Fig. [15] but for run M10J4/5inf at t = 20 = 0.8 Myr. 



